Gellation of rigid filament networks 
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We consider a model for gelation of rigid rods, in which rods that are initially placed at random 
undergo diffusion, and form cross-links when they collide. In the limit of point-like cross-links, 
the number N of croslinks per rod approaches N ~ 3.53. In a model with compliant cross-links of 
maximum length £ c , N(t) increases with time as N(t) oc const. + cL 2 £ c In (i), where c is concentration 
and L is rod length. 

PACS numbers: 61.46.Df,82.70.Kj,87.15.Aa 
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I. INTRODUCTION 

The basic physics of networks of rigid and semiflexible 
filaments has proved to be important to many systems 
of current interest. Carbon nanotubes in a buffered so- 
lution have very strong localized interactions which can 
be modeled as permanent cross-linksi. Similarly, micro- 
tubules and actin filaments in the cytoskeleton have long 
persistence lengths and are cross-linked by a number of 
chemical agents 2 . Still, there is very little knowledge 
concerning the generic properties of cross-linked gels cre- 
ated from long, thin filaments with high bending rigidity. 
There is not even a consensus on the expected connectiv- 
ity of such gels. Without some understanding of connec- 
tivity, it is impossible to formulate ab initio models for 
the electrical and thermal conductivity of these systems, 
much less their rheological properties. 

Previous numerical simulations of the gelation of rigid 
rod networks^ have taken a purely geometrical construc- 
tive approach: In this aproach, straight filaments with 
a nonzero diameter d and a large aspect ratio L/d are 
placed in a unit cell at random, and are assumed to be 
attached only if they overlap. In a three dimensional 
system, the number of cross-links formed in this way de- 
pends directly upon the rod diameter, and vanishes in 
the limit of infinitely thin chains. The situation is, of 
course, very different in a two dimensional model, for 
which this approach can yield a percolation threshhold 
for infinitely thin rods. There is also a rich literature on 
both connectivity and rigidity percolation in networks of 
rods — iSiZiSiS, which rely on such static, constructive ap- 
proaches to building a network. However, when three 
dimensional networks are created in this way, they con- 
tain many near misses, which would be transformed into 
cross-links if the rods were allowed any mobility. The 
number of close approaches just outside the overlap cut- 
off scales with the density of the system. Thus, as we 
will show, this approach vastly underestimates the maxi- 
mal number of connections that can be formed in a three 
dimensional system if the rods are allowed to diffuse. 

Here, we consider a simple dynamical model in which 



gelation instead occurs as the result of Brownian motion 
of thin rods. Rods are initially placed at random and 
then undergo Brownian motion. Pairs of rods are as- 
sumed to irreversibly form cross-links whenever the dis- 
tance of closest approach falls below some "capture ra- 
dius". Throughout this paper, we consider only "flexi- 
ble" cross-links that exert no torque, which thus do not 
introduce a constraint or bias on the angle between cross- 
linked pairs of rods. The simplest variant of such a model 
is one in which rods cross-link when they collide, and 
form a permanent point-like (but rotationally flexible) 
cross-link at the point of collision. In this limit, the cross- 
linking would proceeds to a well defined saturation point: 
The creation of cross-links would continue until the sys- 
tem was mechanically rigid, i.e., until the constraints im- 
posed by the cross-links allowed no further motion of the 
network. 

We may estimate a lower bound on the number of 
cross-links created in this idealized limit of point-like 
cross-links by simple degree of freedom counting: Each 
rod has 5 rigid motion degrees of freedom (three trans- 
lations and two physically relevant rotations, excluding 
axial rotation). Each cross-link introduces three con- 
straints, or removes 3 degrees of freedom, by equating 
the positions of the points along two rods at which those 
rods collide. This argument thus predicts that the to- 
tal ratio of cross-linkers to rods at the rigidity transition 
should be 5/3. Since each cross-link connects two rods, 
the average number of rods to which a randomly chosen 
rod is connected is predicted to be 10/3 = 3.333. 

This estimate is expected to underestimate the num- 
ber of cross-links somewhat, because it assumes that the 
constraints introduced by random cross-links are all lin- 
early independent. It is possible for the system to form 
cross-links that are, in a sense, redundant. The simplest 
example of this occurs when a a rod that is cross-linked 
at two points to a network that would remain rigid if this 
rod was removed. The two cross-links required to immo- 
bilize this rod introduce 6 constraints into the system of 
equations required to describe the system, but remove 
only the 5 degrees of freedom associated with that rod, 
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thus "wasting" one constraint. As a result of this kind 
of redundancy, one might expect the average number of 
attachment points per rod to be slightly higher than 10/3 
(as is found to be the case). 

The idealized limit of point-like cross-links discussed 
above is, however, very difficult to directly simulate, as 
well as being physically unrealistic for many systems of 
interest. The most straightforward way to simulate the 
Brownian motion of a partially gelled system of rods with 
point-like cross-links would be to introduce Langrange 
multiplier forces for each of the cross-links. A Brownian 
dynamics simulation of the system with such cross-links 
would then require the solution of a large system of lin- 
ear equations every time step in order to determine the 
constraint forces exerted by the cross-links. This would 
be prohibitively expensive for all but very small systems 
or short simulations. 

Here, we instead consider a system of slightly compli- 
ant cross-links, in which the distance between the cross- 
link attachment points on a pair of rods is constrained 
to remain less than some maximum cross- link length £ c . 
A network with compliant cross-links retains some flexi- 
bility even when the number of cross-links would arrest 
all motion in a corresponding network of point-like cross- 
links. As a result, the number of cross-links in such a sys- 
tem never truly saturates, but increases more and more 
slowly, as the addition of further cross- links requires more 
and more rare thermal fluctuations. We find that this 
creates a very slow logarithmic growth in the number 
of cross-links with time at late times. We show, how- 
ever, that it is possible to extract information about the 
number of cross-links that would be formed in the ideal- 
ized limit of point-like cross-links by running simulations 
with several values of £ c and extrapolating the results to 

4 = 0. 



II. SIMULATION 

Our simulated systems consist of rigid filaments mod- 
eled by straight line segments of unit length L = 1. The 
initial state of the system is prepared by placing the rods 
in a simulational box with uniformly distributed ran- 
dom positions and orientations. We simulated dimen- 
sionless rod number concentrations c = cL 3 in the range 
2 < c < 2000. Each system is run with at least 1000 
filaments and a box volume > 10L 3 . Periodic boundary 
conditions are enforced. The filaments are moved one at 
a time in random order via an unconstrained Brownian 
dynamics algorithm. For each filament the transverse, 
longitudinal, and rotational displacements were calcu- 
lated separately, with drag coefficients chosen so that the 
transverse and rotation diffusion constants were respec- 
tively 
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FIG. 1: Random distribution of 1000 rigid filaments with den- 
sity c = 10. For filaments which cross the boundaries, both 
periodic images are shown. The diameter of the filaments are 
0.0052 times their length. 



where D\\ is the longitudinal diffusion constant^. 

In the simulations presented here, a permanent cross- 
link is created between two filaments whenever a Brow- 
nian dynamics move would have caused one of the fila- 
ments to cut through the other. Though not presented 
here, we found quantitatively similar results when the 
cross-linking criterion was that two filaments approached 
within a finite capture radius of one another. The cross- 
link connects the two rods at their points of intersection. 
To inhibit bundling, only one cross-link was allowed be- 
tween any two filaments. For dense systems and rigid fil- 
aments, bundling proved not to be an issue, since trapped 
entanglements did not allow filament alignment. 

The constraints of cross-link compliance and topology 
conservation are enforced via a Monte-Carlo acceptance 
criterion, as in RefAi - motions which cause the center- 
lines of neighboring filament to cross or which cause the 
length of any existing cross-link to exceed a maximum 
value are rejected. In some cases we also enforce a hard 
core potential between rods using the same technique. 
Unless otherwise stated, our simulated filaments have 
zero excluded volume. Cross-links imposed no rotational 
constraints about their axes. 

As we show in the next section, we can easily extrap- 
olate from our data to the case of zero radius cross-links 
£ c = 0. We note, however, that direct simulation of cross- 
links that are completely noncompliant to stretching is 
impracticable via any simulational technique. We also 
comment that we are pursuing similar simulations using 
molecular dynamics on bead-chain polymers, but such 
systems are unable to approach the limit of infinite bend 
rigidity. 
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FIG. 2: Time evolution of cross-linking number for £ c — 
0.0013 and, from lowest to highest curves respectively, c = 
100, 150, 200, 300, 400, 600, 800, and 1000. The early time 
evolution of all the plots collapsed when ploted versus c 2 i, 
where time is measured in units of the rotational diffusion 
time l/6D r . 



III. RESULTS 
A. Cross-link number 

Immediately after cross-linking was initiated the num- 
ber of cross-links was found to grow as cVt (see Figure [5]), 
where dimensionless time t is measured in units of the 
rotational diffusion time l/6D r . This power law growth 
ends at a crossover time of order t ~ 100/c 2 . We asso- 
ciate this crossover time with the time necessary for a 
filament to sense its neighbors by diffusion, since the dis- 
tance to nearest neighbors (corresponding to the "cage 
diameter" in the Doi model of entangled rigid rods) is 
proportional to 1 /c. The power law growth was followed 
by a period of slow logarithmic growth. The prefactor of 
the logarithm was observed to increase for larger dimen- 
sionless cross-link radius l c = £ C /L or higher concentra- 
tion c. We simulated many combinations of these param- 
eters in the range 0.0013 < £ c < 0.013 and 5 < c < 2000. 
All simulations were run with dimensionless time-steps 
10 -9 , because we found that the time dependent evolu- 
tion had nearly converged for time-steps on this order 
of magnitude. Obtaining a half decade of logarithmic 
growth on a system of 1000 or more filaments typically 
took 12 days on a single 3.2 GHz Xeon based Linux sys- 
tem. Each data point presented here represents the av- 
erage of four independent systems. 

Figure [3] shows our numerical fits to the scaling data. 
We fit the late time scaling in cross-link number for each 
concentration to the form A\og(c 2 t) + B via parameters 
A and B. Independent fits of A and B to powers of c and 
l c indicated that the overall scaling was consistent with 
a pure dependence on the dimensionless combination c£ c . 
This simple dependence seemed more physically plausi- 
ble, so we restricted our fits to this combination. We 



fit the value of A to a form acl c via parameter a and 
the coefficient B to a form (3 + j3'ct c via the parameters 
(3 and j3' . The parameter (3' may be absorbed into the 
logarithm as a multiplicative constant 7 — exp {(3 1 /a). 

The result of this fitting procedure was the following 
form for the cross-link number: 

N = )3 + ac£ c log(7C 2 f) 
a = 0.0846 ±0.0013 
13 = 3.53 ±0.033 
7 = 4.51 x 10 6 ±2.2 x 10 5 (1) 

We chose to restrict our linear fit of the parameter B to 
the region c£ c < 1.5, which seemed to converge linearly to 
the value B = 3.53 for asymptotically small cl c . Higher 
values measured for c = 200 and 0.078 < 4 < 0.013 
deviated significantly from this linear fit. This deviation 
may result from three body interactions for very large 
values of cross-link radius £ c , or it may simply result from 
our time-step being to large under these very extreme 
conditions. Unfortunately, since these simulations were 
run at the limit of our current computational capability, 
clarification will have to wait for future work. Though 
data is not included here, the cross-link number was seen 
to deviate upward from the fitted linear relation for c < 
10 - we will discuss this discrepency further below. 

In the next section we develop arguments, based on 
additional measurements, for the physical origin of each 
term in Eq. [1] We link the constant term to the rigidifi- 
cation threshold linking number for unstretchable cross- 
links and the term linear in c and l c to localized motion 
within the cross-link length. We attribute the logarith- 
mic growth term to collective motion of the system. The 
logarithm reflects slow, glassy dynamics. 



B. Cross-link statistics 

To shed light on the origins of the cross-linking num- 
ber in Eq [TJ we examine another informative statistic: 
the radial distribution function of unbound neighbors to 
a given rod. It can be shown 12 that for a random spa- 
tial distribution of non- interacting rods, given a test rod, 
the probability that another rod will have a distance of 
closest approach between a radius r and r + dr is given 
by 

P(r)dr = ^cL 2 dr. (2) 

Thus, without cross-linking, the probability of finding a 
neighboring rod with closest approach at any given radius 
is independent of the radius. 

Figure 3] shows the measured distribution for the dis- 
tance of closest approach of unlinked neighbor filaments 
to a given test filament, averaged over all filaments. The 
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FIG. 3: Fitting parameters A and B from a fit to the form j41og (c 2 f) + B for the late time evolution of cross-link number, 
plotted as a function of the product c£ c . Crosses are for runs with constant £ c — 0.0013, stars are for runs with constant 
c = 100, and plus symbols are for runs with constant c = 200. Straight lines are linear fits of A to ac£ c and B to /3 + (5'c£ c 
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FIG. 4: Radial distribution function for distance of closest 
approach of unlinked neighbor filaments. All data is for a 
concentration of cL 3 = 100. The upper and lower data sets 
are respectively for t c /L = 0.0013 and £ C /L = 0.013. Solid 
lines are the predicted distribution from Eq|3]with 6=1. 



measured probability function is well fit by the function 
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where the value of b grows with time. At any given time, 
when P(r) — P c (r) is integrated from zero to infinity it 
gives exactly Eq. [T] for the number of cross-links which 
each filament has captured on average. More spefically, 
the constant term in Eq. [1] comes from integrating the 
exponential term in Eq[31 and all the other terms in Eq.[T] 
comes from itegrating the depleted region between < 
r < bl c in Eq[3] 

We interpret Eq[3]as proof that the extrapolated link- 
age number N = 3.53 is exactly the linkage number re- 
quired for rigidification of line segments in 3-dimensional 



space, which we shall call N r . If the cross-links had no 
stretching compliance (£ c = 0), then after N r cross- links 
per filament had been added the lattice would become 
completely elastically rigid. In this case, there would 
be no further spatial fluctuations, and the cross- linking 
would cease. Prior to the rigidity transition, the filaments 
should diffusively explore their local environment up to 
the point when they have encountered and linked to N r 
other rods on average. Since the probability of finding a 
rod at any given radius is a flat distribution, the proba- 
bility of finding N r rods within radius r is exponential, 
with form 



exp 
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This form correspond exactly to the exponential in Eq[3]if 
we make the identity N r = 3.53. This measured value is 
close to the theoretical lower bound of N r = 10/3 which 
we derived in the introduction. 

We note that as c approaches the value -| x 3.53 ~ 2.2, 
the decay length in Eq. [3] approaches L, so that rod end 
effects will become more important. In this limit the rods 
will interact more like point-like particles. This is con- 
sistent with our observation (not shown) of a deviation 
upward from the logarithmic fit in Eq. [T]for c < 10. 

The effect of finite cross-link length is less straightfor- 
ward to interpret. After the rigidity connectivity N r has 
been reached the only possible fluctuation motion of any 
filament is within the stretch compliance of cross-links to 
its neighbors. The density of neighbors within the cross- 
link length is higher than outside, and without collective 
motion of its neighbors, there is very little room for a 
filament to move around in. Our results suggest that fil- 
aments quickly explore a region of radius r c , producing 
the linear term in Eq[T] At long times, they continue to 
sample and link to neighbors in the space around them, 
increasing their sampling radius at a logarithmic rate. 
Furthermore, the exponential in Eq [2] is shifted by this 
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FIG. 5: Distribution of cross-link separation lengths for late 
time systems with £ c = 0.0013 and c = 150 (lower curve) 
or c = 1000 (upper curve). Straight lines are theoretically 
predicted exponential curve for random cross-link placement, 
disregarding end effects. 
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FIG. 6: Distribution of lengths between filament ends and 
nearest cross-link for late time systems with £ c = 0.0013 and 
c = 150 (lower curve) or c = 1000 (upper curve). 



sampling radius, so that P c {r) goes smoothly to zero at 
the edge of the depletion zone. This is consistent with 
the slow, glassy dynamics of large scale collective motions 
in the presence of deep potential wells. 

Finally, we analyze the statistics of cross-linking. In 
Figure O we plot the distribution of separations between 
neighboring cross-links. As expected, the distribution is 
exponential, with decay constant equal to the average 
density of cross-links per unit length. The exponential 



decay is modified slightly by finite filament length ef- 
fects at separations approaching L. In Figure [6j we plot 
the distribution of lengths of "dangling" ends beyond the 
last cross-link on each filament, which is also exponential. 
These observed distributions are consistent with spatially 
random placements of cross-links along the filament, and 
are qualitatively the same as for flexible polymers. 13 
IV. DISCUSSION 

We have found a new generic number for the maximum 
cross-linking in the limit of vanishing cross-link length, 
iV r = 3.53. We further speculate that this is the as yet 
undiscovered rigidification threshold coordination num- 
ber for a random lattice of rigid, line-like filaments in 3 
dimensions. We have also found that the cross- link num- 
ber depends logarithmically on time and linearly on both 
cross-link length and filament concentration, and we have 
measured the linear coefficient. 

It is worth comparing our results to percolation based 
Monte-Carlo studies without particle dynamics. Foygel 
et al<^ found that the percolation threshold for randomly 
distributed but stationary rods with aspect ratio a oc- 
cured at number density c ~ 0.76a with an average cross- 
linking number N = 1.2. Comparing their simulations to 
ours, their aspect ratio is equivalent to l/(2£ c ). Thus for 
£ c = 0.0013, a typical value in our simulations, connec- 
tivity percolation would not occur below c rs 290, and 
even at this value filaments would have around 1.2 over- 
laps on average. Thus, the addition of dynamics and lo- 
calized bonding between filaments dramatically increases 
the connectivity of the filament suspension. 

An future extension to this work will be the addition 
of a finite filament bend modulus. For large bend modu- 
lus, we expect the transverse filament fluctuations will be 
equivalent to an increase in capture radius. This research 
is under way. It would also be interesting to study the ef- 
fect of the initial configuration on the cross-link number - 
initial alignment of the filaments might greatly increase 
the tendency towards bundling, thereby increasing the 
cross-link number significantly. 



Acknowledgments 

BD thanks Gary Grest and Pieter In't Velt for enlight- 
ening discussions. BD acknowledges partial support from 
the Institute for Mathematics and its Applications with 
funds provided by the National Science Foundation. 



1 L. A. Hough, M. F. Islam, B. Hammouda, A. G. Yodh, P. 
A. Heiney, Nano Letters 6 313 (2006). 

2 B. Alberts, D. Bray, J. Lewis, M. Raff, K. Roberts, and 
J.D. Watson, Molecular Biology of the Cell, 3rd edition 
(Garland, New York, 1994). 



3 I. Balberg, N. Binenbaum, and N. Wagner, Phys. Rev. Lett. 
52, 1465 (1984). 

4 M. Foygel, R. D. Morris, D. Anez, S. French, and V. L. 
Sobolev, Phys. Rev. B 71, 104201 (2005). 

5 D. A. Head, A. J. Levine, and F. C. MacKintosh, Phys. 



6 



Rev. Lett. 91, 108102 (2003). 

6 J. Wilheml and E. Frey, Phys. Rev. Lett. 91, 108103 (2003). 

7 M. Sahaimi and S. Arbabi, Phys. Rev. B 47, 695 (1993), 
M. Sahaimi and S. Arbabi, Phys. Rev. B 47, 703 (1993). 

8 C. Moukarzel and P. M. Duxbury, Phys. Rev. E 59, 2614 
(1999). 

9 D. J. Jacobs and B. Hendrickson, J. Computational 
Physics 137 346 (1997). 

M. Doi and S. F. Edwards, "The Theory of Polymer Dy- 



namics" Oxford University Press, New York (1986). 

11 Shriram Ramanathan and David C. Morse, J. Chera. Phys. 
126, 094906 (2007). 

12 S. Ramanathan, Ph.D. Thesis, University of Minnesota 
(2006). 

13 G. S. Grest and K. Kremer, Macromolecules 23 4994 
(1990). 



